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Abstract 


One-dimensional Navier-Stokes equations are used along with two relaxation 
equations to solve the shock-wave structure for a gas having both rotational and 
vibrational modes of relaxation. The characteristic relaxation times, as well as 
Mach number, are varied in order to note the corresponding nonequilibrium 
effects. Viscosity and heat conduction coefficients are allowed to vary as functions 
of temperature. 

The results presented were obtained by numerical integration of a nonlinear 
set of three differential equations of a boundary value type in which both bound- 
aries exhibited singularities. 
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Shock-Wave Structure of a Gas Having Rotational 
and Vibrational Relaxation 


I. introduction 

The structure of a shock wave is solved by assuming 
that the Navier-Stokes equations and two additional equa- 
tions for rotational and vibrational relaxation rates apply. 
The gas properties can be assumed to be either ideal with 
constant heat capacities and viscosity or real with non- 
constant heat capacities and viscosity as an empirical 
linear function of the appropriate temperature mode. 
This is an extension of the work performed in Ref. 1, but 
differs from their extension, Ref. 2, in that no coupling 
between rotation and vibration occurs. The elimination 
of this coupling mode seems reasonable, since it could 
only be realized through a collision between two mole- 
cules which, in fact, is a translation process. Therefore, 
the following treatment will involve only translational- 
rotational and translational-vibrational coupling. This 
simplification will be useful in future work which will 
include additional nonequihbrium processes such as dis- 
sociation and ionization. It should be noted, however, 
that the validity of the solution is questionable if the 
gradients become too large, since the Navier-Stokes equa- 
tions are not valid for this regime. 


This type of solution to the shock-wave structure is in 
contrast to the more popular technique, the two-step 
model, of assuming a discontinuous shock wave in which 
the translational and rotational temperatures are “shocked” 
and the vibrational temperature remains at the free stream 
condition. A relaxation region then occurs in which the 
vibrational temperature equilibrates with the transla- 
tional and rotational temperature. A shock-structure solu- 
tion of the type reported herein would be useful when 
nonequilibrium effects are known to exist and when the 
density is low, as in the Mars entry problem. 

A number of cases were run which covered an order 
of magnitude in vibrational relaxation time and a factor of 
two in rotational relaxation time. The examples shown 
in this report indicate clearly when the two-step model is 
suflBcient or when the complete shock structure must be 
solved for a solution. An additional advantage of this 
type of solution, aside from accuracy, is that the modified 
Rankine-Hugonoit conditions across the shock wave, 
needed for a solution to a shock-layer problem when the 
shock-wave thickness can not be assumed to be small, 
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are provided. This effect can be seen by imagining that 
the downstream shock conditions for an oblique shock 
wave “shp” along the shock due to the presence of the 
tangential velocity. 

II. Equations of Motion 

The conservation equations of continuity, linear mo- 
mentum, and energy, including viscosity and heat con- 
duction terms, for a one-dimensional steady flow in the 
X direction are 


The continuity equation can be easily integrated, as 

pu — m (4) 

where m is the constant of integration and is equal to the 
mass flux per unit area. 

Using the result of Eq. (4), the linear momentum equa- 
tion can be integrated at once. 

4 du „ 

pu^ + p--^p.-^-P (5) 


d(pu) 

dx 


= 0 


( 1 ) 


dT 

pu c, 


du dp d / 4 du\ _ 
dx dx dx\3^dx) 

dp d / dT\ 4 (du\-^ 

^ dx dx\dx) 3 ^ \dx ) 


( 3 ) 


pU 


/ dT du\ d / dT\ 

dx ^ ^ dxj dx\dx) 


where P is the constant of integration and must be equal 
to Pi + piUi^ or ps + P 2 U 2 since far upstream or far down- 
stream the gardients must vanish and the equation regains 
its Rankine-Hugonoit relationship. 

The integration of the energy equation cannot be per- 
formed directly because of the nonlinearity of the first 
derivative. It can be made integrable, however, if one 
multiplies the hnear momentum equation by u and adds 
it to the energy equation; i.e., 

(£)}=» 


The partitioning of the internal energy of the gas into 
the various modes will be accomplished in Eq. (6). The 
temperature gradient will, therefore, be split into its trans- 
lational, rotational and vibrational parts; i.e.. 


dT 

dTr 

, dTR 

dTy 

^”dx “ 



■ + dx 

dT 

dTv 

dTR 

, dTy 

^ dx 

II 

H . 


dx 


Substituting these into Eq. (6) and then integrating, yield 

( dTr 

pu^ + 2pu {Cpj, Tt + Cj,^ Tr -h Cpy Ty) — 2 < Ay 



where Q is a constant of integration and must be equal 
to the total enthalpy divided by the mass flux per unit 
area. The energy equation, as did the linear momenttun 


equation, has the Rankine-Hugonoit relationship in the 
limi t of far upstream and far downstream from the shock 
wave. 

As can be seen, the number of unknown variables is 
five, whereas there are only three equations. In order to 
obtain two additional equations, it will be assumed that 
the Landau-Teller relaxation model for the rotational and 
vibrational temperatures is valid. 

DTr Tt Tr 

T)t (HR 

for rotation and 

DTv Tt Tv 

Dt (ov 

for vibration, where D/Dt is the material derivative. It 
should be noted that the “bath temperature” is that of 
translation and there is no rotation-vibration coupling 
except indirectly through the translation temperature. 
The characteristic relaxation times for rotation and vibra- 
tion are represented by < 0 ^ and <i>v, respectively. The relaxa- 
tion times can be related to the mean free path L, the 
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mean molecular speed c and the number of collisions N* 
required for the characteristic relaxation time by 


<0i 


NjL 

c 


for rotation, and 

dTy _ Tt - Ty 

dx TtN y [JL 

4 p 


( 12 ) 


The values for oh can be related to viscosity and pressure 
by introducing the following equations: 




fi = pcL, p = pRTt and 
where py. will be taken as equal to Substitution yields 

■7T Ni p 


4 p 


( 10 ) 


Substituting this equation in the relaxation equations, 
and noting that D/Dt = ud/dx for one-dimensional flow, 
we have 


dx ttNb 




4 p 


( 11 ) 


for vibration. Typical values for the number of collisions 
required for equilibration would be of 0(10) for rotation 
and of 0(1000) for vibration. The values of Ny vary 
greatly depending on the type of molecules involved. In 
addition, Ny has a strong pressure dependence and a weak 
temperature dependence. The values for Ny are obtained 
experimentally as no valid theoretical approach has yet 
been devised. A complete discussion of vibrational relaxa- 
tion is presented in Ref. 2. Similar results would probably 
be obtained for rotational relaxation times; but, experi- 
mental values do not exist, except over a limited low 
temperature range. 

The complete set of equations now consist of Eqs. (4, 
5, 9, 11, and 12) and an equation of state. The equation 
of state has been assumed to be that of an ideal gas, 
p = pRTr- The equations will be reduced to a non- 
dimensional form for computational purposes. The non- 
dimensional quantities are introduced as follows: 


mu 




r 


RTr 


p2 ’ 


RTb 

ps 


e = 


m- RTy 
—pi > 



Cy — 


Cvj, 

■r“’ 


CvB 



and 


('Vy 



Also form the follovsdng transport lengths and the non-dimensional counterparts: 


and 




Rr 


mR ’ 


Ly 


Ay 

mR 


JT — 





and 


T*X 



where L* is a reference transport length evaluated at a reference temperature r*. 
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The equations now become 

rv = l (13) 


, 4 / , T* \ dv 

rv^ + 2rv (ctt + + CyS) + 2on 



_ 8 /l^\ ^ _ 
3 \l; )^dy~ 


r" d(j> T T — <jf) 

4 

*r* do T T — 0 

I — „ I -II .1.1 II 

^ dy Lfl "irNy 

4 

<r = Tt 


d<j> T — 

^ dy ~ ttNr 
4 

(16) 

de r-0 

-yS — 

dy irNy 

4 

(17) 

(T = Tt 

(18) 

This set of six equations, of which four are differential 
and two are algebraic, can be reduced to four differential 
equations through the ehmination of a and r by use of 
Eqs. (13) and (18). Carrying this through, one obtains: 

4 dv 

— tV -j— = v^ + t — V 

3 dy 

(19) 

dr d(j> do 

2y.r-^4-2y;,r^4-2yyr-^- 


— v^ + 2v + 2ctt 4- 2cit(j> 4- 2cyd — 

a (20) 


Where Eq. (19) has been used to obtain Eq. (20), Eqs. (16) 
and (17) remain the same since they did not contain o- or t. 


When the viscosity is assumed to be linear^ and in the 
form IX. = (fx.“'/T*) T, where fi* is the viscosity at the refer- 
ence temperature t*, the equations are: 


rv = 1 


(13) 


4 dv 
dy 


— rv^ <7 — 1 


(14) 


dr 

dy 


2 yji T 4- 2yj{T 4" 2yyT 


dcj) 

dy 


dy 


tV 4- rv^ 4- 2rv {cyr 4- 4- Cyd) 


4 " 2<jv 


(15) 


The system of differential equations to be solved are 
improper since the independent variable y does not ap- 
pear anywhere in the equations. It is necessary, therefore, 
to make a change of variable. The convient and logical 
choice is to let v be the independent variable through the 
use of Eq. (19). This will reduce the set of four equations 
in which t, <^, 0, and u are functions of t/ to a set of three 
equations in which t, and 9 are functions of v. Carry- 
ing the transformation through yields 

^ ^ ~ /oi \ 

dv 8 wNk v{v^ + t — v) ^ ' 


^ _ J^_t(t-_0L_ 

dv 3 -rrNy v{v^ + T — V) 


(22) 


and 


dr 2 V^ ( — 4- 2u 4“ 2crT 4- 2cji<j> 4- 2cy<j> — «) — (Syy/TrNy) t (t — 6) — (Sys/TrNs) T (t — (jx) 

dv 3yr « («2 4- T — u) 

This set of differential equations constitutes a two-point boundary value problem with singularities existing at both 
boundaries. 


iThe linear assumption need not be made; however, the formulation would be more difficult. 
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III. Solution of Equations 


The solution of Eqs. (21-23) by numerical means pre- 
sents diflSculties since both the numerator and denomi- 
nator are equal to zero at the initial and final boundary. 
The application of L’Hospital’s rule gives the necessary 
values for 


do 

dv , 


and 


d^ 

dv j 


at the initial boundary. 



dr 

\ _ dv 

1 1 

dv 

3itNr / 
16 

' n .dr 
2oi + 3- 
^ dv 

1 

1^ + Ti 


and 



dr 

1 dv 

1 

dv 

, SrrNy / 

16 

' , dr 

- 1 

1 

^ + Ti 


where the subscript 1 refers to the initial condition. Exam- 
ination of these equations shows that a value of d<f>/dv \ i 
and dO/dv 1 1 can only be found if dr/dv \ i is known. The 
application of L’Kospital’s rule to dr/dv yields a lengthy 
4th degree polynomial in dr/dv | j. The four solutions to 
the polynomial are real, but only one is valid; the other 
three have been extraneously introduced. In all the cases 
tried, the four solutions followed the same pattern repre- 
sented in Fig. 1. Consideration of the physics of the prob- 
lem shows that the value of dr/dv 1 1 has to lay below the 
singularity line and above the dotted straight line joining 
the initial and final value of r. Only one value of dr/dv 1 1 
is between these limits. The sensitivity of the singularity 
makes it necessary to iterate the initial value of dr/dv\i 



Fig. 1 . Initial and final singularities in 
temperature— velocity space 


beyond the accuracy of the value given by the polynomial 
solution. 

The value of dr/dv \i was iterated to 16 significant fig- 
ures, however, the error in r at the final boundary was 
about ±10%. To improve the solution, it was decided to 
take the best two solutions: one falling above and one 
falling below tj, and examine their diEFerences as one 
proceeds away from tj. The At was examined at each 
integration step in order to find the integration step where 
At < e; the next step, however, would make At > €. In 
all cases, e = 10 ® so that At was slightly less than 10"®. 
At this point, r was iterated between the two limits of t 
obtained from the previous “best upper” and ‘Taest lower” 
solutions. This iteration would then yield a new best 
upper and a new best lower solution. This procedure was 
repeated at every point in which r was about to exceed 
10"® until the second boundary was reached. 

Having the solutions of t, <j>, and 0 as a function of v 
through Eqs. (21-23), makes it possible to solve for o- 
and r as a function of v through the use of Eqs. (13) and 
(18). The solutions of t, <j>, 0, a, r, and v may now be 
determined as a function of y through an integration of 
Eq. (19), which was originally used to make the change 
of independent variable from y to v. Since the y = 0 
point is arbitrary, it was taken to be at w = (t>i + v/)/2. 
A discussion of the results of the various cases will be 
given in Section IV. 

IV. Results 

The cases presented are at Mach numbers 1.5, 3.0, and 
10.0 with the number of collisions required for rotational 
equilibration equal to 10 and the number of colhsions 
required for vibrational equilibration equal to 100 and 
500. Two additional cases are presented. The first case 
has a Mach number of 10.0 with the number of collisions 
required for rotational and vibrational equilibration equal 
to 5 and 500, respectively. The second case has a Mach 
number of 3.0 with the number of collisions required for 
rotational and vibrational equilibrium equal to 10 and 
1000, respectively. In all cases, the free-stream conditions 
correspond to those at an altitude of 150,000 ft with 
Ct = 3/2, Cr = 1, Cy = 1, yr = 15/4, ys = 1, and yv = 1, 
and held constant throughout the shock wave. 

The three cases vdth a Mach number of 3.0, Nr = 10, 
and Nv — 100, 500, and 1000 shows the efiEect of vibra- 
tional relaxation time on the shock structure. This is 
shown in Figs. 2, 3, and 4 where the translational, rota- 
tional, and vibrational temperature variation is plotted 
in phase space for Ny = 100, 500, and 1000, respectively. 
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0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

DIMENSIONLESS VELOCITY 


Fig. 2. Temperature vs velocity curves for 
A1 = 3, N« = 10, and Ny = 100 



Fig. 3. Temperature vs velocity curves for 
M=3,Nb= 10, and Ny = 500 


It can be seen that the vibrational temperature rise is 
delayed longer and longer as the number of collisions for 
vibrational equilibration is increased. This effect is more 
noticeable when comparing the Ny = 100 and the 



0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

DIMENSIONLESS VELOCITY 


Fig. 4. Temperature vs velocity curves for 
M = 3, Nr = 10, and Ny = 1000 

Ny = 500 cases as opposed to comparing the Ny = 500 
and the Ny = 1000 cases, since the first comparison has 
a multiplication of 5 in Ny and the second comparison 
has a multiplication of only 2 in Ny. It should also be 
noted that as Ny is increased, the peaks of translational 
and rotational temperature are increased and also delayed 
as is the vibrational temperature rise. Converting from 
the phase space to the physical space by means of Eq. (19) 
results in Figs. 5, 6, and 7. Comparing the corresponding 
figures in phase space and physical space shows the extent 
to which the “relaxation” part of the shock wave is 
expanded in the physical space. 


The density and pressure for M = 3.0, Nr = 10, and 
Ny = 100 in phase space are shown in Fig. 8; Fig. 9 shows 
these parameters when they have been converted to phys- 
ical space. The densities and pressures in phase space 
and physical space, respectively, for M = 3.0, Nr = 10, 
and Nr = 1000 are shown in Figs. 10 and 11. In these 
figures, the effect of converting from the phase space to 
physical space is quite obvious. Again, as in the tempera- 
tures profiles, the relaxation zone has appeared in the 
physical space. One additional feature for density should 
be noted in Fig. 11. The inflection in the density curve 
for Ny = 1000 occurs as the vibrational mode starts its 
change towards equilibrium. This effect is not noticed 
in the density curve for Ny = 100, Fig. 9, since its relax- 
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- 2.0 - 1.0 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 

DIMENSIONLESS DISTANCE 


Fig. 5. Temperature vs distance curves for 
M = 3, Nij = 10, and Ny = 100 



- 2.0 0 2.0 4.0 6.0 8.0 10.0 


DIMENSIONLESS DISTANCE 

Fig. 6. Temperature vs distance curves for 
= 3, Nb = 10, and Ny = 500 


ation starts in the steep part of the profile and is, there- 
fore, not significant. 

The efiFect of Mach number on shock structure can be 
shown in two comparisons. The first comparison is for 
Nk = 10 and Ny = 100 with Mach numbers of 1.5, 3.0, 
and 10.0. Temperature profiles in phase space are shown 
in Figs. 12, 2, and 13, respectively. It is interesting to note 
that the effect of increasing only the Mach number is 
similar to that of increasing only the relaxation time for 
vibrational equilibration. These same temperature pro- 
files, when they are presented in physical space, are shown 
in Figs. 14, 5, and 15. The increased overshoot of trans- 
lational temperature, as the Mach number is increased, 
is the interesting feature; this is due to the increased 
amount of kinetic energy that must be absorbed in the 
early part of the shock by the translational temperature. 
The density and pressure profiles in phase space are shown 
in Figs. 16, 8, and 17; Figs. 18, 9, and 19 show these pro- 
files in physical space. These curves show no secondary 
relaxation region as did the case for Ny — 1000, with 
M = 3.0, and N® = 10 as shown in Fig. 11 (special note 
should be taken of the density profile). The second com- 
parison is for the same value of Nr and the same Mach 
number values with Ny now equal to 500. The tempera- 
ture profiles in phase space for these cases are shown in 
Figs. 20, 3, and 21. The same effect, as was pointed out 
for the cases with Ny = 100, is also present, although not 



DIMENSIONLESS DISTANCE 


Fig. 7. Temperature vs distance curves for 
M=3,Nr= 10, and Ny = 1000 

to such a great extent. The temperature profiles in physi- 
cal space are shown in Figs. 22, 6, and 23. Density and 
pressure profiles are shown in phase space in Figs. 24 
and 25 (the profiles for the M = 3.0 case have been 
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DIMENSIONLESS PRESSURE AND DENSITY DIMENSIONLESS PRESSURE AND DENSITY 


5.0 



Fig. 8. Pressure and density vs velocity curves 
for M=3,Nk = 10, and Nr = 100 



»ui„. ,1 , , .J I — I 

0.2 0.4 0.6 0.8 1.0 


DIMENSIONLESS VELOCITY 

Fig. 10. Pressure and density vs velocity curves 
for M = 3, Nb = 10, and Ny = 1000 



Fig. 9. Pressure and density vs distance curves 
for M = 3, Nr = 10, and Nr ~ 100 


omitted) and in physical space in Figs. 26 and 27. In con- 
trast to the cases in the first comparison, a relaxation 
region has appeared and can be seen in both the density 
and pressure profiles in physical space. 


One additional case in which Nr = 5, Ny = 500, and 
a Mach number of 10.0 is presented in Figs. 28 and 29 for 
temperature profiles in phase space and physical space, 
respectively. Comparing these illustrations with Figs. 21 
and 23 for Ny = 500 and Mach number of 10.0 but with 
Nr = 10 shows a marked diflFerence in the shock wave 
thickness. This is shown in the physical space where for 
Nr — 5 the rotational temperature equilibrates with trans- 
lational temperature at t/ ^ 0.82; whereas for Nr = 10, the 
equilibration does not take place until y ^ 1.30. The 
equilibration for vibrational temperature is correspond- 
ingly longer for IVr = 10 as compared to Nr = 5. 


V. Conclusions 

The cases presented show that for all relaxation times 
in the vibrational mode, there was a significant rise in 
the vibrational temperature before the rotational tem- 
perature equilibrated with the translational temperature. 
As can be seen in Figs. 5, 6, and 7 for Nr = 100, 500, and 
1000 the percentage rise in vibrational temperature before 
rotation-translation equilibration was 45, 16.5, and 13.5%, 
respectively. For this reason, a shock-structure analysis 
is preferred to the “two-step” model when accuracy is 
required to better than say 10%. This last statement; of 
course, only holds for Ny ^500. For Ny < 500, the two- 
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5.0 



DIMENSIONLESS DISTANCE 

Fig. 1 1 . Pressure and density vs distance curves 
for A1 = 3, Nb = 10, and Ny = 1000 

step model is completely unsatisfactory and the only 
recourse is a shock-structure analysis. 

The values of Nb were not varied much since they are 
generally accepted to be in the neighborhood of from 
5 to 10 collisions. For a given value of Ny, the two-step 
model will be better as Nr is decreased and will, con- 
versely, be worse as Nr is increased. In this respect, Nr 
plays a large role as to whether or not the two-step model 
is applicable to a given problem. 

As stated earlier, care should be taken when using the 
shock-structure analysis for high Mach numbers and high 



0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 

DIMENSIONLESS VELOCITY 


Fig. 12. Temperature vs velocity curves for 
M = 1 .5, Nb = 10, and Ny = 100 

densities, since large gradients could introduce errors into 
the results through the Navier-Stokes equations. Whether 
significant errors occur is a function of the values for Nr 
and Ny as well as the Mach number; therefore, each case 
must be examined individually. 

It is clear that this problem could be extended to in- 
clude dissociation, and possibly ionization, by the use of 
an augmented equation of state along with a rate equa- 
tion for dissociation, or ionization. For example, Lighthill’s 
ideal dissociating gas could be used for the equation of 
state in the dissociation region and a similar equation 
could be applied in the region in which ionization occurs. 
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DIMENSIONLESS TEMPERATURE DIMENSIONLESS TEMPERATURE 



Fig. 13. Temperature vs velocity curves for 
Ai = 10, Nr = 10, and Nv = 100 



- 10.0 0 10.0 20.0 30.0 


DIMENSIONLESS DISTANCE 

Fig. 14. Temperature vs distance curves for 
M = 1.5, Nr = 10, and Ny = 100 



Fig. 15. Temperature vs distance curves for 
M = 10, Nr = 10, and Ny = 100 



0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 


DIMENSIONLESS VELOCITY 

Fig. 16. Pressure and density vs velocity curves 
for M = 1.5, Nr = 10, and Ny = 100 
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DIMENSIONLESS PRESSURE AND DENSITY DIMENSIONLESS PRESSURE AND DENSITY 



DIMENSIONLESS VELOCITY 


Fig. 17. Pressure and density vs velocity curves 
for M = 10, Nb = 10, and Ny = 100 



- 10.0 - 5.0 0 5.0 10.0 15.0 20.0 25.0 30.0 

DIMENSIONLESS DISTANCE 


Fig. 1 8. Pressure and density vs distance curves 
for M = 1.5, = 10, and Ny = 100 






DIMENSIONLESS TEMPERATURE DIMENSIONLESS TEMPERATURE 



0.1 0.2 0.3 0.4 0.5 0.5 0.7 0.8 0.9 1.0 


DIMENSIONLESS VELOCITY 

Fig. 21. Temperature vs velocity curves for 
M = 10, Nb = 10, and Ny = 500 
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Fig. 23. Temperature vs distance curves for 
M = 10, Nb = 10, and Ny = 500 
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Fig. 22. Temperature vs distance curves for 
A1 = 1.5, Nb = 10, and Ny = 500 
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Fig. 24. Pressure and density vs velocity curves 
for M = 1 .5, Nb = 10, and Ny = 500 
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DIMENSIONLESS PRESSURE AND DENSITY DIMENSIONLESS PRESSURE AND DENSITY 



Fig. 25. Pressure and density vs velocity curves 
for A1 = 10, Nb = 10, and Ny = 500 



Fig. 26. Pressure and density vs distance curves 
for M = 1.5, Nb = 10, and Ny = 500 
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Fig. 29. Temperature vs distance curves for 
M = 10, N„ = 5, and Ny = 500 
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